rm(list=ls())
suppressPackageStartupMessages({
library(rjags)
library(bayesplot)
library(jagsUI)
library(glue)
library(loo)
library(knitr)
library(kableExtra)
}
)
Setting up environment, importing the dataset and visualising a graph of the two time series we will have to fit.
rm(list=ls())
# If set to TRUE higher order models will also be trained
intense = TRUE
# New Family Houses Sold: United States
# Source: https://fred.stlouisfed.org/series/HSN1F
data <- read.csv("gdp_inflation.csv",header=T)
data$DATE <- as.Date(data$DATE)
# Plotting time series for visualisation purposes
plot(
data$DATE[1:305],
data$GDP_PC1[1:305],
type="l",
xlab="",
ylab="",
main="Percent change YoY (quarterly)"
)
lines(
data$DATE[1:305],
data$CPIAUCSL_PC1[1:305],
type="l",
col="red"
)
legend(x="bottom",legend=c("GDP","CPI"),col=c("black","red"),lty=1)
At this point, we setup our environment by converting the values in both time series into numericals, and then splitting them into training and testing. Finally, we define some useful variables which we will be using through out our script.
# Dataset extraction
x_axis <- data$DATE[1:305]
series_GDP <- as.numeric(data$GDP_PC1[1:305])
series_CPI <- as.numeric(data$CPIAUCSL_PC1[1:305])
# Parameter setup (RMK: we assume both time series of same length)
TRAIN_PERC <- 0.15
N_tot <- length(series_GDP)
N_test <- floor(TRAIN_PERC * N_tot)
N_train <- N_tot - N_test
To further inspect our data we can also plot the respective auto-correlation functions:
acf(series_GDP, lag.max=100)
acf(series_CPI, lag.max=100)
The auto-correlation functions of the two proceses do not clearly match any given noticeable pattern, hence a more sophisticated model selection will be needed, where multiple models will be tested and evaluated one against the other.
For our convenience, we define a few utility functions which will be useful later. The first one automatises the interactions with JAGS, the second one will draw plots of both the whole time-series and its out-of-sample portion only, the third one will specifically plot a time-series along with the anomalies which have been detected through out its history, and so on.
fit_JAGS <- function(
model_string, # JAGS model string
data_list, # List of input values for JAGS
to_save, # Vector with parameter names which should be saved
n_iter=5000, # Number of iterations. Should be greater than n_burnin
n_adapt=1000, # JAGS n.adapt internal parameter
n_chain=1, # Number of chains to run
n_burnin=1000, # Burn-in value
n_thin=1, # Thinning value
verbose=FALSE # Displays JAGS logging / output (turn on for debugging)
) {
output <- jags(
model.file=textConnection(model_string),
data=data_list,
parameters.to.save=to_save,
n.adapt=n_adapt,
n.iter=n_iter,
n.chains=n_chain,
n.burnin=n_burnin,
n.thin=n_thin,
verbose=verbose
)
return(output)
}
visualise <- function(
x_axis, # Values to be displayed on the x-axis of the plots
original_series, # Original time series with all samples
output, # Output from JAGS, with in-sample and out-of-sample predictions
title="", # Title to display on the graphs
keep=5 # Training sampled to keep in the test graph
) {
# Length of total series
N <- length(x_axis)
# Plotting entire graph (training-set + predictions)
# In-sample predicted values with their 2.5%-97.5% quantiles
yp_in_mean <- output$mean$yp_in
yp_in_q1 <- output$q2.5$yp_in
yp_in_q2 <- output$q97.5$yp_in
in_length <- length(yp_in_mean)
# Out-of-sample predicted values with their 2.5%-97.5% quantiles and 25%-75% quantiles
yp_out_mean <- output$mean$yp_out
yp_out_q1 <- output$q2.5$yp_out
yp_out_q2 <- output$q97.5$yp_out
yp_out_q3 <- output$q25$yp_out
yp_out_q4 <- output$q75$yp_out
# Removing potentially non-finite values from quantiles
yp_in_q1[is.infinite(yp_in_q1)] <- NA
yp_in_q2[is.infinite(yp_in_q2)] <- NA
yp_out_q1[is.infinite(yp_out_q1)] <- NA
yp_out_q2[is.infinite(yp_out_q2)] <- NA
yp_out_q3[is.infinite(yp_out_q3)] <- NA
yp_out_q4[is.infinite(yp_out_q4)] <- NA
# Computing extreme values for first graph
min_y <- min(original_series, yp_in_q1, yp_out_q1)
max_y <- max(original_series, yp_in_q2, yp_out_q2)
# Plotting REAL time series (as points only), both training set and test-set
plot(
x_axis,
original_series,
main=title, col="red",
xlab="time", ylab="price",
ylim=c(min_y, max_y)
)
lines(
x_axis,
original_series,
col="red"
)
# Plotting a vertical line signalling end of training set and horizontal
# Showing the computer mean
abline(v=x_axis[in_length], col="magenta", lwd=1, lty=2)
abline(h=output$mean$m0, col="green", lwd=1, lty=2)
# Plotting in-sample predictions as well as 2.5%-97.5% quantiles
points(x_axis[1:in_length], yp_in_mean, col="blue", pch="*")
lines(x_axis[1:in_length], yp_in_q1, col="blue", lwd=1.5)
lines(x_axis[1:in_length], yp_in_q2, col="blue", lwd=1.5)
# Plotting out-of-sample predicted values as well as 2.5%-97.5%
# quantiles and 25%-75% quantiles
points(x_axis[(in_length+1):N], yp_out_mean, col="orange", pch="*")
lines(x_axis[(in_length+1):N], yp_out_q1, col="orange", lwd=1.5)
lines(x_axis[(in_length+1):N], yp_out_q2, col="orange", lwd=1.5)
lines(x_axis[(in_length+1):N], yp_out_q3, col="orange", lwd=1.5, lty=2)
lines(x_axis[(in_length+1):N], yp_out_q4, col="orange", lwd=1.5, lty=2)
# Plotting onyl out-of-sample portion of graph (with last 5 training samples)
# Re-computing extreme values for second graph
min_y <- min(original_series[(in_length+1 - keep):N], yp_out_q1)
max_y <- max(original_series[(in_length+1 - keep):N], yp_out_q2)
# Plotting REAL time series (as points only), both training set and test-set
plot(
x_axis[(in_length+1 - keep):N],
original_series[(in_length+1 - keep):N],
main=title, col="red",
xlab="time", ylab="price",
ylim=c(min_y, max_y)
)
lines(
x_axis[(in_length+1-keep):N],
original_series[(in_length+1-keep):N],
col="red"
)
# Plotting a vertical line signalling end of training set
abline(v=x_axis[in_length+1], col="magenta", lwd=1, lty=2)
abline(h=output$mean$m0, col="green", lwd=1, lty=2)
# Plotting out-of-sample predicted values as well as 2.5%-97.5%
# quantiles and 25%-75% quantiles
points(x_axis[(in_length+1):N], yp_out_mean, col="orange", pch="*")
lines(x_axis[(in_length+1):N], yp_out_q1, col="orange", lwd=1.5)
lines(x_axis[(in_length+1):N], yp_out_q2, col="orange", lwd=1.5)
lines(x_axis[(in_length+1):N], yp_out_q3, col="orange", lwd=1.5, lty=2)
lines(x_axis[(in_length+1):N], yp_out_q4, col="orange", lwd=1.5, lty=2)
}
# Function in charge of visualising VAR models (requires custom rendering).
# This acts as a wrapper function to the standard visualise(...) function, defined above.
visualise_VAR = function(
x_axis, # Values to be displayed on the x-axis of the plots
original_series_gdp, # Original GDP time series with all samples
original_series_infl, # Original INFLATION time series with all samples
output, # Output from VAR JAGS, with in-sample and out-of-sample predictions
title="VAR", # Title to display on the graphs
keep=5 # Training sampled to keep in the test graph
) {
# Marginalise GDP
output_gdp = list()
output_gdp$mean$yp_in = output$mean$yp_in_gdp
output_gdp$q2.5$yp_in = output$q2.5$yp_in_gdp
output_gdp$q97.5$yp_in = output$q97.5$yp_in_gdp
output_gdp$mean$yp_out = output$mean$yp_out_gdp
output_gdp$q2.5$yp_out = output$q2.5$yp_out_gdp
output_gdp$q97.5$yp_out = output$q97.5$yp_out_gdp
output_gdp$q25$yp_out = output$q25$yp_out_gdp
output_gdp$q75$yp_out = output$q75$yp_out_gdp
output_gdp$mean$m0 = output$mean$m_gdp
# Marginalise CPI
output_infl = list()
output_infl$mean$yp_in = output$mean$yp_in_infl
output_infl$q2.5$yp_in = output$q2.5$yp_in_infl
output_infl$q97.5$yp_in = output$q97.5$yp_in_infl
output_infl$mean$yp_out = output$mean$yp_out_infl
output_infl$q2.5$yp_out = output$q2.5$yp_out_infl
output_infl$q97.5$yp_out = output$q97.5$yp_out_infl
output_infl$q25$yp_out = output$q25$yp_out_infl
output_infl$q75$yp_out = output$q75$yp_out_infl
output_infl$mean$m0 = output$mean$m_infl
# Call visualisation on marginalisations
visualise(
x_axis,
original_series_gdp,
output_gdp,
title=glue("GDP series, {title}"),
keep
)
visualise(
x_axis,
original_series_infl,
output_infl,
title=glue("CPI series, {title}"),
keep
)
}
# Function in charge of visualising the anomaly points on top of its corresponding time series
visualise_anomalies <- function(
x_axis, # X-axis to use for the plot
original_series, # Original time series to display
high_var_indices, # List of indices corresponding to high variance points
med_var_indices, # List of indices corresponding to medium variance points
title="", # Name of the time series used
start_point=1 # Value after which plotting should start (in order to obtain close up graphs)
) {
# Plot actual time series
N <- length(x_axis)
plot(
x_axis[start:N],
original_series[start:N],
main=glue("Anomalies for {title}"),
t="l", xlab="time", ylab="price"
)
# Plotting strong anomaly points
points(x_axis[high_var_indices], original_series[high_var_indices], col="red")
for (high_anomaly in high_var_indices) {
abline(v=x_axis[high_anomaly], col="red", lwd=1.5)
}
# Plotting medium anomaly points
points(x_axis[med_var_indices], original_series[med_var_indices], col="orange")
for (medium_anomaly in med_var_indices) {
abline(v=x_axis[medium_anomaly], col="orange", lty=3, lwd=1.5)
}
}
# Function in charge of classifying anomaly points given their delta value after running the anomaly detection simulation
classify_anomalies <- function(
output, # JAGS simulation output
N_tot, # Total number of samples
title="", # Name of time series to display
high_value_thr=0.2, # High variance threshold value used
med_value_thr=0.60 # Medium variance threshold value used
) {
# Obtaining precision, variance and delta
tau <- output$mean$tau
sigma2 <- 1 / tau
delta <- output$mean$d
# Classifying the values of d
high_var_indices <- seq(1:(N_tot-1))[delta < high_value_thr]
med_var_indices <- seq(1:(N_tot-1))[high_value_thr < delta & delta < med_value_thr]
low_var_indices <- seq(1:(N_tot-1))[delta > med_value_thr]
# Logging counts
n_high <- length(high_var_indices)
n_med <- length(med_var_indices)
print(
glue("Found {n_high} strong and {n_med} medium outliers out {N_tot} total samples for {title}")
)
# Returning the indices
return(list(
high_var_indices = high_var_indices,
med_var_indices = med_var_indices,
low_var_indices = low_var_indices
))
}
# Function in charge of plotting the anomaly point delta values with their relative threshold lines
display_classification <- function(
x_axis, # X-axis to use for the plot
delta, # Sequence of delta values
high_var_indices, # List of indices corresponding to high variance points
med_var_indices, # List of indices corresponding to medium variance points
low_var_indices, # List of indices corresponding to low variance points
title="", # Name of the time series used
high_value_thr=0.2, # High variance threshold value used
med_value_thr=0.60 # Medium variance threshold value used
) {
# Plotting the posterior values of delta and classifying them
plot(
x_axis[low_var_indices], delta[low_var_indices],
xlab="time",
ylab="Posterior delta values",
main=glue("{title} anomaly samples"),
xlim=c(min(x_axis), max(x_axis)),
ylim=c(0, 1)
)
points(x_axis[high_var_indices], delta[high_var_indices], col="red")
points(x_axis[med_var_indices], delta[med_var_indices] , col="orange")
abline(h=high_value_thr, col="red")
abline(h=med_value_thr, col="orange")
}
# K is the number of effective params
compute_metrics <- function(output, series, N_train, k) {
############################################################
###################### IN-SAMPLE ###########################
############################################################
# In-sample: DIC, WAIC, BIC
ll_mat <- output$sims.list$loglik
DIC_val <- compute_dic_marginal(ll_mat)
waic_res <- suppressWarnings( waic(ll_mat) )
WAIC_val <- waic_res$estimates["waic","Estimate"]
#BIC computation
#Sum of the log-likelihoods for each draw
ll_sum <- rowSums(ll_mat)
# maximum of the sum: estimate of log p(y|theta_hat), that is, log-lik_max
lhat <- max(ll_sum)
# BIC_val <- -2 * lhat + k * log(N_train)
BIC_val <- -2 * lhat + k * log(N_train - output$mean$train_metrics_skip)
# BIC Approx
# BIC_val <- DIC_val + k * (log(N_train) - 2)
############################################################
#################### OUT-OF-SAMPLE #########################
############################################################
y_pred <- if (is.null(output$mean$yp_onestep_out)) output$mean$yp_out else output$mean$yp_onestep_out
y_true <- series[(N_train+1):length(series)]
e <- y_pred - y_true
# MSE
MSE_val <- mean((e)^2)
#MAE
MAE_val <- mean(abs(e))
MAPE_val <- mean(abs(e / y_true)) * 100
#SMAPE
SMAPE_val <- mean(2 * abs(e) / (abs(y_true) + abs(y_pred))) * 100
#MASE
scale <- mean(abs(diff(series[1:N_train])))
MASE_val <- MAE_val / scale
list(
DIC = round(DIC_val, 1),
WAIC = round(WAIC_val, 1),
BIC = round(BIC_val, 1),
MSE = round(MSE_val, 3),
MAE = round(MAE_val, 3),
MAPE = round(MAPE_val, 1),
SMAPE = round(SMAPE_val, 1),
MASE = round(MASE_val, 3)
)
}
compute_dic_marginal <- function(ll_mat) {
dev <- -2 * ll_mat
dev_sum <- rowSums(dev)
Dbar <- mean(dev_sum)
pD <- 0.5 * var(dev_sum)
DIC <- Dbar + pD
return(DIC)
}
# Appends metrics row to provided metrics_df and returns it
# If no arguments are passed, generates and returns a new empty metrics_df
metrics_df_append <- function(metrics_df = NA, new_mets = NA, name="...") {
if (anyNA(metrics_df)) {
return(data.frame(
Model = character(),
DIC = numeric(),
WAIC = numeric(),
BIC = numeric(),
MSE = numeric(),
MAE = numeric(),
MAPE = numeric(),
SMAPE = numeric(),
MASE = numeric(),
stringsAsFactors = FALSE
))
} else if (!anyNA(metrics_df) && !all(is.na(new_mets))) {
return(rbind(
metrics_df,
data.frame(
Model = name,
DIC = new_mets$DIC,
WAIC = new_mets$WAIC,
BIC = new_mets$BIC,
MSE = new_mets$MSE,
MAE = new_mets$MAE,
MAPE = new_mets$MAPE,
SMAPE = new_mets$SMAPE,
MASE = new_mets$MASE,
stringsAsFactors = FALSE
))
)
}
return(metrics_df)
}
At this point, we define all our JAGS model strings so as to be able
to reuse them multiple times through out the script. We start with with
defining the model for a generic \(\text{AR}(n)\) sequence. The
auto-regressive order \(n\) is passed
to the model by means of the order variable.
\[ y(t) = \left[ \sum_{i=1}^{n} \alpha_i y(t-i) \right] + \eta(t) \qquad \eta \sim \mathcal{N}(0, \sigma^2) \]
For any order \(n ≥ 1\) the following model will correctly implement an \(\text{AR}(n)\) process, which allows us to automatically test multiple AR processes easily and at once.
# AR(n) string
modelAR_string <- "
model {
############################################################
######################## LIKELIHOOD ########################
############################################################
# 1. Bootstrapping the AR process
mu[1:order] <- y[1:order]
yp_in[1:order] <- y[1:order]
# 2. Recursive step of the AR process
for (t in (order+1):N_train) {
mu[t] <- inprod(alpha_rev, y[(t-order):(t-1)]) + m0
y[t] ~ dnorm(mu[t], tau)
yp_in[t] ~ dnorm(mu[t], tau)
}
for (t in (train_metrics_skip+1):N_train) {
loglik[t-train_metrics_skip] <- logdensity.norm(y[t], mu[t], tau)
}
############################################################
######################## PREDICTION ########################
############################################################
# Init helper array z
z[1:order] <- y[(N_train-order+1):N_train]
for (t in (order+1):(order+N_test)) {
mu_z[t] <- inprod(alpha_rev, z[(t-order):(t-1)]) + m0
z[t] ~ dnorm(mu_z[t], tau)
}
yp_out <- z[(order+1):(order+N_test)]
############################################################
############ ONE STEP AHEAD PREDICTION #####################
############################################################
for (t in 1:N_test) {
mu_onestep[t] <- inprod(alpha_rev, y[(N_train+t-order):(N_train+t-1)]) + m0
yp_onestep_out[t] ~ dnorm(mu_onestep[t], tau)
}
############################################################
########################## PRIOR ###########################
############################################################
m0 ~ dnorm(0.0, 1.0E-4)
tau ~ dgamma(0.1, 10)
for (i in 1:order) {
alpha[i] ~ dunif(-1.5, 1.5)
alpha_rev[order - i + 1] <- alpha[i]
}
# Defining converstion from precision to variance here
sigma2 <- 1 / tau
}"
We now define a generic \(\text{MA}(n)\) process. The memory order
\(n\) is passed to the model by means
of the order variable.
\[ y(t) = \left[ \sum_{i=1}^n \beta_i \eta(t-i) \right] + \eta(t) \qquad \eta \sim \mathcal{N}(0, \sigma^2) \]
For any order \(n ≥ 1\) the following model will correctly implement a \(\text{MA}(n)\) process, which allows us to automatically test multiple MA processes easily and at once.
modelMA_string <- "
model {
############################################################
######################## LIKELIHOOD ########################
############################################################
e[1:order] = rep(0, order)
yp_in[1:order] <- y[1:order]
# 2. Recursive step of the MA process
for (t in (order+1):N_train) {
mu[t] <- m0 + inprod(beta_rev, e[(t-order):(t-1)])
y[t] ~ dnorm(mu[t], tau)
yp_in[t] ~ dnorm(mu[t], tau)
e[t] <- y[t] - mu[t]
}
for (t in (train_metrics_skip+1):N_train) {
loglik[t-train_metrics_skip] <- logdensity.norm(y[t], mu[t], tau)
}
############################################################
######################## PREDICTION ########################
############################################################
for (t in 1:N_test) {
e[N_train+t] ~ dnorm(0, tau)
mu_out[t] <- m0 + inprod(beta_rev, e[(N_train+t-order):(N_train+t-1)])
yp_out[t] ~ dnorm(mu_out[t], tau)
}
############################################################
############ ONE STEP AHEAD PREDICTION #####################
############################################################
e_onestep[1:order] <- e[(N_train-order+1):N_train]
for (t in 1:N_test) {
mu_onestep[t] <- m0 + inprod(beta_rev, e_onestep[t:(t+order-1)])
yp_onestep_out[t] ~ dnorm(mu_onestep[t], tau)
e_onestep[t+order] <- y[N_train+t] - mu_onestep[t]
}
############################################################
########################## PRIOR ###########################
############################################################
m0 ~ dnorm(0.0, 1.0E-4)
tau ~ dgamma(0.01, 0.01)
for (i in 1:order){
beta[i] ~ dnorm(0, 1/10)
beta_rev[i] <- beta[order-i+1]
}
# Defining converstion from precision to variance here
sigma2 <- 1 / tau
}"
Finally, the following model will merge the two previous models defining an \(\text{ARMA}(p, q)\) model.
\[ y(t) = \left[ \sum_{i=1}^p \alpha_i y(t-i) \right] + \left[ \sum_{i=1}^q \beta_i \eta(t-i) \right] + \eta(t) \qquad \eta \sim \mathcal{N}(0, \sigma^2). \]
The parameters order_ar and order_ma allow
us to specify the two orders when simulating the model so as to be able
to test all possible combinations of auto-regressive and moving-average
processes.
modelARMA_string <- "
model {
############################################################
######################## LIKELIHOOD ########################
############################################################
e[1:order] <- rep(0, order)
mu[1:order] <- y[1:order]
yp_in[1:order] <- y[1:order]
# 2. Recursive step of the MA process
for (t in (order+1):N_train) {
mu_ar[t] <- inprod(alpha_rev, y[(t-order):(t-1)])
mu_ma[t] <- inprod(beta_rev, e[(t-order):(t-1)])
mu[t] <- m0 + mu_ar[t] + mu_ma[t]
y[t] ~ dnorm(mu[t], tau)
yp_in[t] ~ dnorm(mu[t], tau)
e[t] <- y[t] - mu[t]
}
for (t in (train_metrics_skip+1):N_train) {
loglik[t-train_metrics_skip] <- logdensity.norm(y[t], mu[t], tau)
}
############################################################
######################## PREDICTION ########################
############################################################
z[1:order] <- y[(N_train-order+1):N_train]
e_z[1:order] <- e[(N_train-order+1):N_train]
for (t in (order+1):(order+N_test)) {
e_z[t] ~ dnorm(0, tau)
mu_ar_out[t] <- inprod(alpha_rev, z[(t-order):(t-1)])
mu_ma_out[t] <- inprod(beta_rev, e_z[(t-order):(t-1)])
mu_z[t] <- m0 + mu_ar_out[t] + mu_ma_out[t]
z[t] ~ dnorm(mu_z[t], tau)
}
yp_out <- z[(order+1):(order+N_test)]
############################################################
############ ONE STEP AHEAD PREDICTION #####################
############################################################
e_onestep[1:order] <- e[(N_train-order+1):N_train]
for (t in 1:N_test) {
mu_ar_onestep[t] <- inprod(alpha_rev, y[(N_train+t-order):(N_train+t-1)])
mu_ma_onestep[t] <- inprod(beta_rev, e_onestep[t:(t+order-1)])
mu_onestep[t] <- m0 + mu_ar_onestep[t] + mu_ma_onestep[t]
yp_onestep_out[t] ~ dnorm(mu_onestep[t], tau)
e_onestep[t+order] <- y[N_train+t] - mu_onestep[t]
}
############################################################
########################## PRIOR ###########################
############################################################
m0 ~ dnorm(0.0, 1.0E-4)
tau ~ dgamma(0.01, 0.01)
for (i in 1:order_ar) {
alpha[i] ~ dunif(-1.5, 1.5)
}
for (i in (order_ar+1):order) {
alpha[i] <- 0
}
for (i in 1:order) {
alpha_rev[i] <- alpha[order-i+1]
}
for (i in 1:order_ma) {
beta[i] ~ dnorm(0, 1/10)
}
for (i in (order_ma+1):order) {
beta[i] <- 0
}
for (i in 1:order) {
beta_rev[i] <- beta[order-i+1]
}
# Defining converstion from precision to variance here
sigma2 <- 1 / tau
}"
At this point we can start modelling our two time series using the stochastic processes for which we defined the previous model strings. We will start with the GDP time-series and then focus on the CPI time-series. Finally, we will also attempt to model a correlated process using a \(\text{VAR}(n)\) model, where both time-series will be modelled together at once.
# Metrics dataframe setup
metrics_df_GDP <- metrics_df_append()
We will, at first, start off with some \(\text{AR}(n)\) models, with increasing values of \(n\). We will use the obtained results ot decide whether or not the selected model may apply to the problem at hand, and we will also evaluate which regressive order is best.
to_save <- c(
"alpha",
"m0",
"sigma2",
"yp_in",
"yp_out",
"yp_onestep_out",
"train_metrics_skip",
"loglik"
)
# Trying our AR(n) from n=1 to max_order
max_order <- if (intense) 9 else 2
for (order in 1:max_order) {
data_list <- list(
"y"=series_GDP,
"order"=order,
"N_train"=N_train,
"N_test"=N_test,
"train_metrics_skip"=50
)
# Perform simulations and visualise output
output <- fit_JAGS(
modelAR_string,
data_list,
to_save
)
visualise(x_axis, series_GDP, output, title=glue("GDP series, AR({order})"))
# Update metrics
mets <- compute_metrics(output, series_GDP, N_train, k=order+2)
metrics_df_GDP <- metrics_df_append(metrics_df_GDP, mets, glue("AR({order})"))
}
We will now also fit some \(\text{MA}(n)\) processes (for now \(1 \leq n \leq 5\)) to evaluate whether anyone might be better.
to_save <- c(
"beta",
"m0",
"sigma2",
"yp_in",
"yp_out",
"yp_onestep_out",
"train_metrics_skip",
"loglik"
)
# Trying our MA(n) from n=1 to max_order
max_order <- if (intense) 5 else 2
for (order in 1:max_order) {
data_list <- list(
"y"=series_GDP,
"order"=order,
"N_train"=N_train,
"N_test"=N_test,
"train_metrics_skip"=50
)
# Perform simulations and visualise output
output <- fit_JAGS(
modelMA_string,
data_list,
to_save
)
visualise(x_axis, series_GDP, output, title=glue("GDP series, MA({order})"))
# Updating metrics
mets <- compute_metrics(output, series_GDP, N_train, k=order+2)
metrics_df_GDP <- metrics_df_append(metrics_df_GDP, mets, glue("MA({order})"))
}
Finally, we combine the auto-regressive and moving-average proccesses and try to model our series using an \(\text{ARMA}(p, q)\) process. For now, we will attempt modelling for all combinations of \(1 \leq p \leq 4\) and \(1 \leq q \leq 4\).
to_save <- c(
"alpha",
"beta",
"m0",
"sigma2",
"yp_in",
"yp_out",
"yp_onestep_out",
"train_metrics_skip",
"loglik"
)
# Trying our AR(n) from n=1 to max_order
max_order_ar <- if (intense) 4 else 2
max_order_ma <- if (intense) 4 else 2
for (order_ma in 1:max_order_ma) {
for (order_ar in 1:max_order_ar) {
data_list <- list(
"y"=series_GDP,
"N_train"=N_train,
"N_test"=N_test,
"order"=max(order_ar, order_ma),
"order_ar"=order_ar,
"order_ma"=order_ma,
"train_metrics_skip"=50
)
# Perform simulations and visualise output
output <- fit_JAGS(
modelARMA_string,
data_list,
to_save
)
visualise(x_axis, series_GDP, output, title=glue("GDP series, ARMA({order_ar}, {order_ma})"))
# Updating metrics
mets <- compute_metrics(output, series_GDP, N_train, k=order_ar+order_ma+2)
metrics_df_GDP <- metrics_df_append(metrics_df_GDP, mets, glue("ARMA({order_ar}, {order_ma})"))
}
}
# Metrics dataframe setup
metrics_df_CPI <- metrics_df_append()
We continue our modelling moving to the CPI time-series. We will proceede as for the GDP, trying our a few \(\text{AR}(n)\) models, then moving to some \(\text{AR}(n)\) and finally, putting both together, some \(\text{ARMA}(p, q)\) models.
to_save <- c(
"alpha",
"m0",
"sigma2",
"yp_in",
"yp_out",
"yp_onestep_out",
"train_metrics_skip",
"loglik"
)
# Trying our AR(n) from n=1 to max_order
max_order <- if (intense) 9 else 2
for (order in 1:max_order) {
data_list <- list(
"y"=series_CPI,
"order"=order,
"N_train"=N_train,
"N_test"=N_test,
"train_metrics_skip"=50
)
# Perform simulations and visualise output
output <- fit_JAGS(
modelAR_string,
data_list,
to_save
)
visualise(x_axis, series_CPI, output, title=glue("CPI series, AR({order})"))
if (order == 9) {
output_ar9_cpi = output
}
# Update metrics
mets <- compute_metrics(output, series_CPI, N_train, k=order+2)
metrics_df_CPI <- metrics_df_append(metrics_df_CPI, mets, glue("AR(",order,")"))
}
We will now also fit some \({MA}(n)\) processes (for now \(1 \leq n \leq 5\)) to evaluate whether anyone might be better.
to_save <- c(
"beta",
"m0",
"sigma2",
"yp_in",
"yp_out",
"yp_onestep_out",
"train_metrics_skip",
"loglik"
)
# Trying our MA(n) from n=1 to max_order
max_order <- if (intense) 5 else 2
for (order in 1:max_order) {
data_list <- list(
"y"=series_CPI,
"order"=order,
"N_train"=N_train,
"N_test"=N_test,
"train_metrics_skip"=50
)
# Perform simulations and visualise output
output <- fit_JAGS(
modelMA_string,
data_list,
to_save
)
visualise(x_axis, series_CPI, output, title=glue("CPI series, MA({order})"))
# Updating metrics
mets <- compute_metrics(output, series_CPI, N_train, k=order+2)
metrics_df_CPI <- metrics_df_append(metrics_df_CPI, mets, glue("MA(",order,")"))
}
Finally, we combine the auto-regressive and moving-average proccesses and try to model our series using an \(\text{ARMA}(p, q)\) process. For now, we will attempt modelling for all combinations of \(1 \leq p \leq 4\) and \(1 \leq q \leq 4\).
to_save <- c(
"alpha",
"beta",
"m0",
"sigma2",
"yp_in",
"yp_out",
"yp_onestep_out",
"train_metrics_skip",
"loglik"
)
# Trying our AR(n) from n=1 to max_order
max_order_ar <- if (intense) 4 else 2
max_order_ma <- if (intense) 4 else 2
for (order_ma in 1:max_order_ma) {
for (order_ar in 1:max_order_ar) {
data_list <- list(
"y"=series_CPI,
"N_train"=N_train,
"N_test"=N_test,
"order"=max(order_ar, order_ma),
"order_ar"=order_ar,
"order_ma"=order_ma,
"train_metrics_skip"=50
)
# Perform simulations and visualise output
output <- fit_JAGS(
modelARMA_string,
data_list,
to_save
)
visualise(x_axis, series_CPI, output, title=glue("CPI series, ARMA({order_ar}, {order_ma})"))
# Updating metrics
mets <- compute_metrics(output, series_CPI, N_train, k=order_ar+order_ma+2)
metrics_df_CPI <- metrics_df_append(metrics_df_CPI, mets, glue("ARMA({order_ar}, {order_ma})"))
}
}
Finally, we will try to fit generic order \(\text{VAR}(n)\) models, which have the following form:
\[ \vec{y}(t) = \vec{\mu}_0 + \left[ \sum_{i=1}^n A_i \vec{y}(t-i) \right] + \vec{\eta}(t) \qquad \forall i \; A_i \in \mathbb{M}(2 \times 2, \mathbb{R}) \qquad \vec{\eta} \sim \mathcal{N}(\vec{0}, \Sigma) \]
Where \(\Sigma\) is the covariance
matrix of our white noise (now an \(\mathbb{R}^2\) vector) and \(\vec{y}(t) \in \mathbb{R}^2\) now contains
both the GDP and CPI time-series, which will be modelled together in a
correlated fashion. By passing the parameter order during
simulation we can specify the auto-regressive order \(n\) of the multivariate process.
The idea behind this model is to attempt to model not only the two time-series independently, but also in a correlated way, as we work under the assumption that the two may influence each other.
# Prior fragment for Wishart distributions (positive definite covariance)
prior_wishart <- "
# Wishart prior for a full covariance matrix
R[1,1] <- 1
R[2,2] <- 1
R[1,2] <- 0
R[2,1] <- 0
tau ~ dwish(R, 3) # df=3 is the minimum for a 2x2 matrix
"
# Prior fragment for independent Gamma distributions (diagonal covariance)
prior_diagonal <- "
# Independent Gamma priors for diagonal elements of precision matrix
tau[1,1] ~ dgamma(0.01, 0.01)
tau[2,2] ~ dgamma(0.01, 0.01)
tau[1,2] <- 0
tau[2,1] <- 0
"
modelVAR_string = "
model {
############################################################
######################## LIKELIHOOD ########################
############################################################
mu[1:2,1:order] = y[,1:order]
yp_in[1:2,1:order] = y[,1:order]
for (t in (order+1):N_train) {
mu_tmp[1:2, t, 1] <- m0
for (j in 1:order) {
mu_tmp[1:2,t,j+1] <- mu_tmp[1:2, t, j] + A[,,j] %*% y[,(t-j)]
}
mu[1:2,t] <- mu_tmp[1:2,t,order+1]
y[1:2,t] ~ dmnorm(mu[1:2,t], tau)
yp_in[1:2,t] ~ dmnorm(mu[1:2,t], tau)
}
for (t in (train_metrics_skip+1):N_train) {
loglik_gdp[t-train_metrics_skip] <- logdensity.norm(y[1,t], mu[1,t], 1 / sigma2[1,1])
loglik_infl[t-train_metrics_skip] <- logdensity.norm(y[2,t], mu[2,t], 1 / sigma2[2,2])
}
############################################################
######################## PREDICTION ########################
############################################################
z[1:2,1:order] = y[,(N_train-order+1):N_train]
for (t in (order+1):(order+N_test)) {
mu_z[1:2, t, 1] <- m0
for (j in 1:order) {
mu_z[1:2,t,j+1] <- mu_z[1:2, t, j] + A[,,j] %*% z[,(t-j)]
}
z[1:2,t] ~ dmnorm(mu_z[1:2,t,order+1], tau)
}
yp_out[1:2,1:N_test] = z[1:2,(order+1):(order+N_test)]
############################################################
############ ONE STEP AHEAD PREDICTION #####################
############################################################
for (t in 1:N_test) {
mu_onestep_tmp[1:2, t, 1] <- m0
for (j in 1:order) {
mu_onestep_tmp[1:2,t,j+1] <- mu_onestep_tmp[1:2, t, j] + A[,,j] %*% y[,(N_train+t-j)]
}
mu_onestep[1:2,t] <- mu_onestep_tmp[1:2,t,order+1]
yp_onestep_out[1:2,t] ~ dmnorm(mu_onestep[1:2,t], tau)
}
############################################################
########################## PRIOR ###########################
############################################################
for (i in 1:2) {
for (j in 1:2) {
for (o in 1:order) {
A[i, j, o] ~ dunif(-1.2, 1.2)
# A[i, j, o] ~ dnorm(0, 10)
}
}
}
# Dynamically inserted prior for tau
{{prior_for_tau}}
# __if UNCORRELATED_ERRORS__
#
# tau[1,1] ~ dgamma(0.01, 0.01)
# tau[2,2] ~ dgamma(0.01, 0.01)
# tau[1,2] <- 0
# tau[2,1] <- 0
# __else__
#
# tau ~ dwish(R, 3) # dwish(R=ScaleMatrix... Identity, df=degrees of freedom)
m_gdp ~ dnorm(0.0, 1.0E-4)
m_infl ~ dnorm(0.0, 1.0E-4)
m0 <- c(m_gdp, m_infl)
sigma2 <- inverse(tau)
############################################################
######################## EXTRACTION ########################
############################################################
# Expose them for output in R
loglik_gdp_out <- loglik_gdp
loglik_infl_out <- loglik_infl
yp_in_gdp <- yp_in[1,]
yp_in_infl <- yp_in[2,]
yp_out_gdp <- yp_out[1,]
yp_out_infl <- yp_out[2,]
yp_onestep_out_gdp <- yp_onestep_out[1,]
yp_onestep_out_infl <- yp_onestep_out[2,]
}"
# Setting up parameters
y_var <- matrix(c(series_GDP, series_CPI), nrow=2, byrow=TRUE)
order <- 1
use_wishart <- TRUE
if (use_wishart) {
prior_for_tau = prior_wishart
} else {
prior_for_tau = prior_diagonal
}
final_modelVAR_string <- glue(modelVAR_string, .open="{{", .close="}}")
# Uncomment for debugging inspection
#print(glue(modelVAR_string, .open="{{", .close="}}"))
# Defining JAGS parameters
to_save <- c(
"A",
"m_gdp",
"m_infl",
"sigma2",
"yp_in_gdp",
"yp_in_infl",
"yp_out_gdp",
"yp_out_infl",
"yp_onestep_out_gdp",
"yp_onestep_out_infl",
"train_metrics_skip",
"loglik_gdp", # marginal log-likelihood for GDP
"loglik_infl" # marginal log-likelihood for CPI
)
data_list <- list(
"y"=y_var,
"order"=order,
"N_train" = N_train,
"N_test" = N_test,
"train_metrics_skip"=50
)
# Fitting VAR model and visualising it
jags_output_VAR = fit_JAGS(
final_modelVAR_string,
data_list,
to_save,
n_iter=6000,
n_adapt=1000,
n_chain=1,
n_burnin=2000
)
visualise_VAR(
x_axis,
series_GDP,
series_CPI,
jags_output_VAR,
title=glue("VAR({order})"),
keep=5
)
# Wrapper for GDP-only metrics
output_gdp <- list(
mean = list(
yp_out = jags_output_VAR$mean$yp_out_gdp,
yp_onestep_out = jags_output_VAR$mean$yp_onestep_out_gdp,
train_metrics_skip=jags_output_VAR$mean$train_metrics_skip
),
sims.list = list(loglik = jags_output_VAR$sims.list$loglik_gdp)
)
# Wrapper for CPI-only metrics
output_infl <- list(
mean = list(
yp_out = jags_output_VAR$mean$yp_out_infl,
yp_onestep_out = jags_output_VAR$mean$yp_onestep_out_infl,
train_metrics_skip=jags_output_VAR$mean$train_metrics_skip
),
sims.list = list(loglik = jags_output_VAR$sims.list$loglik_infl)
)
# Compute univariate metrics for GDP and CPI as well
# RMK: just for visualisation purposes: can't compare marginalised metrics with univariate metrics
mets_VAR_gdp <- compute_metrics(output_gdp, series_GDP, N_train, k = 4*order + 2 + 3)
mets_VAR_infl <- compute_metrics(output_infl, series_CPI, N_train, k = 4*order + 2 + 3)
metrics_df_GDP <- metrics_df_append(metrics_df_GDP, mets_VAR_gdp, glue("VAR_univar_GDP({order})"))
metrics_df_CPI <- metrics_df_append(metrics_df_CPI, mets_VAR_infl, glue("VAR_univar_GDP({order})"))
At this point, we display all so far computer metrics for all models and both time-series, so as to be able to perform model selection objectively.
# GDP table
kable(
metrics_df_GDP,
caption = "Univariate model metrics — GDP series",
digits = c(1,1,1,1,3,3,1,1,3),
format = "html",
table.attr = 'cellpadding="10" cellspacing="0"'
) %>%
kable_styling(full_width = F, position = "center") %>%
row_spec(0, background = "#000000", color = "#FFFFFF") %>% # intestazione
row_spec(1:nrow(metrics_df_GDP), background = "#000000", color = "#FFFFFF") # righe dati
| Model | DIC | WAIC | BIC | MSE | MAE | MAPE | SMAPE | MASE |
|---|---|---|---|---|---|---|---|---|
| AR(1) | 744.7 | 705.6 | 694.4 | 7.971 | 1.378 | 75.0 | 26.5 | 1.195 |
| AR(2) | 690.0 | 667.2 | 666.5 | 12.174 | 1.828 | 144.0 | 31.9 | 1.585 |
| AR(3) | 690.5 | 666.3 | 666.8 | 11.280 | 1.731 | 134.9 | 30.5 | 1.501 |
| AR(4) | 682.5 | 663.8 | 673.2 | 11.885 | 1.756 | 135.7 | 30.0 | 1.523 |
| AR(5) | 644.2 | 630.5 | 647.0 | 8.451 | 1.543 | 124.5 | 27.0 | 1.338 |
| AR(6) | 647.2 | 632.9 | 653.0 | 8.151 | 1.522 | 121.9 | 26.9 | 1.320 |
| AR(7) | 654.9 | 636.7 | 658.0 | 8.193 | 1.532 | 121.0 | 25.9 | 1.329 |
| AR(8) | 660.6 | 637.8 | 660.1 | 8.373 | 1.590 | 125.6 | 28.3 | 1.379 |
| AR(9) | 645.4 | 626.7 | 655.6 | 7.683 | 1.494 | 124.1 | 26.6 | 1.296 |
| MA(1) | 846.8 | 835.5 | 838.8 | 9.079 | 2.077 | 114.2 | 42.6 | 1.801 |
| MA(2) | 792.0 | 780.5 | 786.8 | 15.539 | 2.488 | 130.8 | 45.0 | 2.157 |
| MA(3) | 659.1 | 638.7 | 643.1 | 5.097 | 1.316 | 100.7 | 27.8 | 1.142 |
| MA(4) | 620.8 | 604.1 | 613.0 | 5.912 | 1.344 | 140.0 | 27.6 | 1.166 |
| MA(5) | 652.9 | 620.5 | 627.9 | 7.048 | 1.388 | 133.0 | 25.5 | 1.203 |
| ARMA(1, 1) | 697.9 | 673.7 | 672.0 | 11.378 | 1.880 | 136.7 | 34.7 | 1.630 |
| ARMA(2, 1) | 684.6 | 664.0 | 669.9 | 11.637 | 1.789 | 138.8 | 31.4 | 1.551 |
| ARMA(3, 1) | 675.6 | 654.3 | 662.1 | 11.163 | 1.799 | 165.1 | 32.8 | 1.560 |
| ARMA(4, 1) | 690.3 | 661.6 | 667.4 | 9.795 | 1.701 | 141.0 | 31.9 | 1.475 |
| ARMA(1, 2) | 683.2 | 671.7 | 661.2 | 6.573 | 1.457 | 115.2 | 28.5 | 1.263 |
| ARMA(2, 2) | 646.4 | 631.6 | 639.4 | 7.314 | 1.559 | 136.2 | 30.5 | 1.352 |
| ARMA(3, 2) | 632.2 | 613.6 | 626.9 | 7.601 | 1.567 | 132.5 | 30.3 | 1.359 |
| ARMA(4, 2) | 649.3 | 624.2 | 638.6 | 7.902 | 1.560 | 129.1 | 30.4 | 1.353 |
| ARMA(1, 3) | 592.2 | 585.0 | 598.6 | 6.245 | 1.295 | 142.4 | 24.9 | 1.123 |
| ARMA(2, 3) | 598.2 | 571.4 | 582.2 | 5.323 | 1.099 | 130.7 | 22.1 | 0.953 |
| ARMA(3, 3) | 594.0 | 579.7 | 598.5 | 5.837 | 1.204 | 137.2 | 23.5 | 1.044 |
| ARMA(4, 3) | 611.2 | 590.1 | 610.5 | 5.958 | 1.251 | 136.6 | 24.2 | 1.085 |
| ARMA(1, 4) | 601.2 | 582.0 | 595.5 | 6.111 | 1.260 | 138.7 | 24.2 | 1.093 |
| ARMA(2, 4) | 614.7 | 583.4 | 593.7 | 6.259 | 1.285 | 142.0 | 24.8 | 1.115 |
| ARMA(3, 4) | 608.7 | 587.2 | 606.7 | 6.126 | 1.267 | 138.2 | 24.4 | 1.098 |
| ARMA(4, 4) | 611.3 | 589.3 | 615.6 | 5.957 | 1.253 | 136.9 | 24.3 | 1.086 |
| VAR_univar_GDP(1) | 740.9 | 704.1 | 725.0 | 8.013 | 1.381 | 76.1 | 26.9 | 1.198 |
# CPI table
kable(
metrics_df_CPI,
caption = "Univariate model metrics — CPI series",
digits = c(1,1,1,1,3,3,1,1,3),
format = "html",
table.attr = 'cellpadding="10" cellspacing="0"'
) %>%
kable_styling(full_width = F, position = "center") %>%
row_spec(0, background = "#000000", color = "#FFFFFF") %>%
row_spec(1:nrow(metrics_df_CPI), background = "#000000", color = "#FFFFFF")
| Model | DIC | WAIC | BIC | MSE | MAE | MAPE | SMAPE | MASE |
|---|---|---|---|---|---|---|---|---|
| AR(1) | 569.2 | 552.4 | 550.4 | 0.689 | 0.607 | 202.3 | 33.5 | 0.899 |
| AR(2) | 562.1 | 548.3 | 552.0 | 0.580 | 0.592 | 217.8 | 35.0 | 0.877 |
| AR(3) | 557.6 | 546.6 | 554.1 | 0.559 | 0.578 | 200.9 | 37.8 | 0.857 |
| AR(4) | 542.0 | 532.6 | 546.4 | 0.566 | 0.571 | 169.0 | 36.4 | 0.847 |
| AR(5) | 509.7 | 502.2 | 520.3 | 0.508 | 0.575 | 152.1 | 36.4 | 0.853 |
| AR(6) | 506.9 | 498.2 | 518.2 | 0.511 | 0.580 | 183.8 | 36.5 | 0.860 |
| AR(7) | 510.5 | 499.8 | 520.9 | 0.508 | 0.577 | 191.6 | 36.3 | 0.855 |
| AR(8) | 511.2 | 500.4 | 526.3 | 0.515 | 0.581 | 190.2 | 36.0 | 0.861 |
| AR(9) | 482.6 | 468.2 | 499.0 | 0.447 | 0.528 | 263.3 | 34.3 | 0.783 |
| MA(1) | 800.3 | 792.5 | 797.2 | 1.887 | 1.181 | 913.3 | 52.7 | 1.752 |
| MA(2) | 737.9 | 725.9 | 727.1 | 1.321 | 0.972 | 651.8 | 47.1 | 1.441 |
| MA(3) | 571.7 | 560.5 | 571.4 | 0.589 | 0.594 | 441.3 | 34.0 | 0.880 |
| MA(4) | 526.1 | 507.8 | 515.0 | 0.573 | 0.579 | 423.0 | 35.5 | 0.859 |
| MA(5) | 551.3 | 526.7 | 536.1 | 0.467 | 0.531 | 336.0 | 34.9 | 0.787 |
| ARMA(1, 1) | 556.3 | 546.8 | 551.4 | 0.606 | 0.605 | 217.9 | 35.9 | 0.898 |
| ARMA(2, 1) | 546.5 | 527.7 | 527.5 | 0.619 | 0.632 | 276.1 | 41.3 | 0.937 |
| ARMA(3, 1) | 546.7 | 540.8 | 556.6 | 0.561 | 0.577 | 201.1 | 37.3 | 0.855 |
| ARMA(4, 1) | 511.5 | 499.6 | 509.1 | 0.537 | 0.609 | 268.6 | 40.4 | 0.903 |
| ARMA(1, 2) | 541.8 | 532.4 | 534.1 | 0.529 | 0.524 | 143.7 | 34.5 | 0.777 |
| ARMA(2, 2) | 543.3 | 526.2 | 531.4 | 0.638 | 0.639 | 304.3 | 41.6 | 0.948 |
| ARMA(3, 2) | 501.7 | 494.5 | 509.2 | 0.473 | 0.508 | 141.1 | 35.0 | 0.753 |
| ARMA(4, 2) | 504.5 | 497.5 | 514.1 | 0.566 | 0.605 | 291.2 | 40.3 | 0.898 |
| ARMA(1, 3) | 477.8 | 468.0 | 479.2 | 0.424 | 0.534 | 304.3 | 35.0 | 0.792 |
| ARMA(2, 3) | 467.7 | 456.5 | 471.7 | 0.401 | 0.520 | 259.5 | 34.8 | 0.770 |
| ARMA(3, 3) | 436.9 | 421.3 | 435.7 | 0.378 | 0.496 | 230.3 | 35.0 | 0.735 |
| ARMA(4, 3) | 467.2 | 447.7 | 464.2 | 0.436 | 0.525 | 266.5 | 36.0 | 0.779 |
| ARMA(1, 4) | 464.1 | 449.1 | 462.4 | 0.405 | 0.493 | 222.1 | 33.5 | 0.731 |
| ARMA(2, 4) | 478.8 | 449.0 | 458.2 | 0.426 | 0.484 | 221.6 | 33.2 | 0.717 |
| ARMA(3, 4) | 476.6 | 452.2 | 462.1 | 0.410 | 0.489 | 219.6 | 33.5 | 0.725 |
| ARMA(4, 4) | 457.2 | 440.5 | 464.1 | 0.449 | 0.525 | 214.4 | 36.4 | 0.779 |
| VAR_univar_GDP(1) | 550.5 | 537.6 | 567.1 | 0.721 | 0.636 | 200.4 | 39.3 | 0.944 |
NOTICE: The MAPE metric exceeds 100% because the series (quarterly or monthly GDP growth rates) has true values (denominators) that are very small or close to zero. This results in an amplified percentage.
We will now display some trace plots for some parameters of some of the fitted models and relative posterior density plots to make sure the MCMC samples are not getting stuck in some local regions, and are behaving well overall.
n_chain = if (intense) 3 else 1
to_save <- c(
"alpha",
"m0",
"sigma2"
)
data_list <- list(
"y"=series_GDP,
"order"=9,
"N_train"=N_train,
"N_test"=N_test,
"train_metrics_skip"=50
)
output <- fit_JAGS(
modelAR_string,
data_list,
to_save,
n_thin=100,
n_adapt=2000,
n_chain=n_chain,
n_iter=100000,
verbose=TRUE
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 251
## Unobserved stochastic nodes: 352
## Total graph size: 1901
##
## Initializing model
##
## Adaptive phase, 2000 iterations x 3 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 1000 iterations x 3 chains
##
##
## Sampling from joint posterior, 99000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.
color_scheme_set("teal")
mcmc_trace(output$samples[,1:11])
mcmc_dens(output$samples[,1:11])
mcmc_intervals(output$samples[,1:11])
n_chain = if (intense) 3 else 1
to_save <- c(
"beta",
"m0",
"sigma2"
)
data_list <- list(
"y"=series_GDP,
"order"=3,
"N_train"=N_train,
"N_test"=N_test,
"train_metrics_skip"=50
)
output <- fit_JAGS(
modelMA_string,
data_list,
to_save,
n_thin=100,
n_adapt=2000,
n_chain=n_chain,
n_iter=100000,
verbose=TRUE
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 257
## Unobserved stochastic nodes: 397
## Total graph size: 2265
##
## Initializing model
##
## Adaptive phase, 2000 iterations x 3 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 1000 iterations x 3 chains
##
##
## Sampling from joint posterior, 99000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.
mcmc_trace(output$samples[,1:5])
mcmc_dens(output$samples[,1:5])
mcmc_intervals(output$samples[,1:5])
## --- 1. Support Function: matrix [p × q] ---------------------------
make_matrix <- function(df, metric, P, Q){
M <- matrix(NA_real_, nrow = P, ncol = Q) # row = p, column = q
for(i in seq_len(nrow(df))){
nums <- as.integer(regmatches(df$Model[i],
gregexpr("\\d+", df$Model[i]))[[1]])
if(length(nums) == 2){
p <- nums[1]; q <- nums[2]
M[p, q] <- df[i, metric]
}
}
M
}
## --- 2. function for plotting axis -------------------------
plot_heat <- function(mat, title){
image(x = 1:ncol(mat),
y = 1:nrow(mat),
z = mat,
col = colorRampPalette(
c("navy","deepskyblue","yellow","orange","red"))(100),
axes = FALSE, xlab = "", ylab = "", main = title, useRaster = TRUE)
contour(1:ncol(mat), 1:nrow(mat), mat,
add = TRUE, lwd = 0.7, drawlabels = FALSE)
axis(1, at = 1:ncol(mat), labels = 1:ncol(mat), cex.axis = 0.75) # q
axis(2, at = 1:nrow(mat), labels = 1:nrow(mat), las = 1,
cex.axis = 0.75) # p
box()
}
## --- 3. parameters & palette --------------------------------------------
max_p <- max_order_ar # max order AR
max_q <- max_order_ma # max order MA
metrics <- c("DIC","WAIC","BIC","MSE",
"MAE","MAPE","SMAPE","MASE")
## --- 4. heat-map 2×4 for axis p↔q --------------------------------------
op <- par(mfrow = c(2,4), mar = c(2.2,2.2,2.4,0.6))
for(met in metrics){
mat <- make_matrix(metrics_df_GDP, met, max_p, max_q)
plot_heat(mat, paste("ARMA –", met))
}
par(op)
n_chain = if (intense) 3 else 1
to_save <- c(
"alpha",
"beta",
"m0",
"sigma2"
)
order_ar <- 2
order_ma <- 3
data_list <- list(
"y"=series_GDP,
"N_train"=N_train,
"N_test"=N_test,
"order"=max(order_ar, order_ma),
"order_ar"=order_ar,
"order_ma"=order_ma,
"train_metrics_skip"=50
)
# Perform simulations and visualise output
output <- fit_JAGS(
modelARMA_string,
data_list,
to_save,
n_thin=50,
n_adapt=2000,
n_chain=n_chain,
n_iter=50000,
verbose=TRUE
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 257
## Unobserved stochastic nodes: 399
## Total graph size: 2964
##
## Initializing model
##
## Adaptive phase, 2000 iterations x 3 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 1000 iterations x 3 chains
##
##
## Sampling from joint posterior, 49000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.
color_scheme_set("teal")
mcmc_trace(output$samples[,1:8])
mcmc_dens(output$samples[,1:8])
mcmc_intervals(output$samples[,1:8])
n_chain = if (intense) 3 else 1
to_save <- c(
"alpha",
"m0",
"sigma2"
)
data_list <- list(
"y"=series_CPI,
"order"=9,
"N_train"=N_train,
"N_test"=N_test,
"train_metrics_skip"=50
)
output <- fit_JAGS(
modelAR_string,
data_list,
to_save,
n_thin=100,
n_adapt=2000,
n_chain=n_chain,
n_iter=100000,
verbose=TRUE
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 251
## Unobserved stochastic nodes: 352
## Total graph size: 1901
##
## Initializing model
##
## Adaptive phase, 2000 iterations x 3 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 1000 iterations x 3 chains
##
##
## Sampling from joint posterior, 99000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.
color_scheme_set("teal")
mcmc_trace(output$samples[,1:11])
mcmc_dens(output$samples[,1:11])
mcmc_intervals(output$samples[,1:11])
n_chain = if (intense) 3 else 1
to_save <- c(
"beta",
"m0",
"sigma2"
)
data_list <- list(
"y"=series_CPI,
"order"=4,
"N_train"=N_train,
"N_test"=N_test,
"train_metrics_skip"=50
)
output <- fit_JAGS(
modelMA_string,
data_list,
to_save,
n_thin=100,
n_adapt=2000,
n_chain=n_chain,
n_iter=100000,
verbose=TRUE
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 256
## Unobserved stochastic nodes: 397
## Total graph size: 2261
##
## Initializing model
##
## Adaptive phase, 2000 iterations x 3 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 1000 iterations x 3 chains
##
##
## Sampling from joint posterior, 99000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.
mcmc_trace(output$samples[,1:6])
mcmc_dens(output$samples[,1:6])
mcmc_intervals(output$samples[,1:6])
## --- 1. Support Function: matrix [p × q] ---------------------------
make_matrix <- function(df, metric, P, Q){
M <- matrix(NA_real_, nrow = P, ncol = Q)
for(i in seq_len(nrow(df))){
nums <- as.integer(regmatches(df$Model[i],
gregexpr("\\d+", df$Model[i]))[[1]])
if(length(nums) == 2){
p <- nums[1]; q <- nums[2]
M[p, q] <- df[i, metric]
}
}
M
}
## --- 2. function for plotting axis -------------------------
plot_heat <- function(mat, title){
image(x = 1:ncol(mat),
y = 1:nrow(mat),
z = mat,
col = colorRampPalette(
c("navy","deepskyblue","yellow","orange","red"))(100),
axes = FALSE, xlab = "", ylab = "", main = title, useRaster = TRUE)
contour(1:ncol(mat), 1:nrow(mat), mat,
add = TRUE, lwd = 0.7, drawlabels = FALSE)
axis(1, at = 1:ncol(mat), labels = 1:ncol(mat), cex.axis = 0.75) # q
axis(2, at = 1:nrow(mat), labels = 1:nrow(mat), las = 1,
cex.axis = 0.75) # p
box()
}
## --- 3. parameters & palette --------------------------------------------
max_p <- max_order_ar # max order AR
max_q <- max_order_ma # max order MA
metrics <- c("DIC","WAIC","BIC","MSE",
"MAE","MAPE","SMAPE","MASE")
## --- 4. heat-map 2×4 for axis p↔q --------------------------------------
op <- par(mfrow = c(2,4), mar = c(2.2,2.2,2.4,0.6))
for(met in metrics){
mat <- make_matrix(metrics_df_CPI, met, max_p, max_q)
plot_heat(mat, paste("ARMA –", met))
}
par(op)
n_chain = if (intense) 3 else 1
to_save <- c(
"alpha",
"beta",
"m0",
"sigma2"
)
order_ar <- 3
order_ma <- 3
data_list <- list(
"y"=series_CPI,
"N_train"=N_train,
"N_test"=N_test,
"order"=max(order_ar, order_ma),
"order_ar"=order_ar,
"order_ma"=order_ma,
"train_metrics_skip"=50
)
# Perform simulations and visualise output
output <- fit_JAGS(
modelARMA_string,
data_list,
to_save,
n_thin=100,
n_adapt=2000,
n_chain=n_chain,
n_iter=100000,
verbose=TRUE
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 257
## Unobserved stochastic nodes: 400
## Total graph size: 2965
##
## Initializing model
##
## Adaptive phase, 2000 iterations x 3 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 1000 iterations x 3 chains
##
##
## Sampling from joint posterior, 99000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.
color_scheme_set("teal")
mcmc_trace(output$samples[,1:8])
mcmc_dens(output$samples[,1:8])
mcmc_intervals(output$samples[,1:8])
to_save <- c(
"A",
"m0",
"sigma2"
)
data_list <- list(
"y"=y_var,
"order"=1,
"N_train"=N_train,
"N_test"=N_test,
"train_metrics_skip"=50
)
output <- fit_JAGS(
final_modelVAR_string,
data_list,
to_save,
n_thin=3,
n_adapt=2000,
n_chain=n_chain,
n_iter=6000,
verbose=TRUE
)
##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 259
## Unobserved stochastic nodes: 356
## Total graph size: 2735
##
## Initializing model
##
## Adaptive phase, 2000 iterations x 3 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 1000 iterations x 3 chains
##
##
## Sampling from joint posterior, 5000 iterations x 3 chains
##
##
## Calculating statistics.......
##
## Done.
mcmc_trace(output$samples[,1:10])
mcmc_dens(output$samples[,1:10])
mcmc_intervals(output$samples[,1:10])
Based on the calculated metrics, the models that seem to perform best are \(\text{ARMA}(2, 3)\) and \(\text{ARMA}(3, 3)\), both for the GDP and inflation time series.
An honorable mention must be given to the \(\text{AR}(9)\) model for the CPI time series, where the out-of-sample prediction for the next year and a half, given only the in-sample data, is almost perfectly aligned with the actual values, as shown in the following plot.
if (intense) {
visualise(x_axis, series_CPI, output_ar9_cpi, title=glue("CPI series, AR(9)"))
}
In this last section we will examine anomalies and outliers which are present in our two time series. So far, we performed model selection and predictions of future values for the series, however we could clearly see that both the training and test set presented samples which deviated quite heavily from their neighboring values. This can be cause of disturbances during the learning phase and also a reduced accuracy during forecasting.
The objective of this section is to detect said outliers and visualise them so as to clearly show which samples may be problematic during fitting and out-of-sample prediction.
In order to do this, for each timestamp \(t\), we consider the incremental values of our series, each of which is to be treated as an independent random variable:
\[ \Delta_t := y_{t} - y_{t-1} \stackrel{iid}{\sim} \mathcal{N}(\mu, \sigma_t^2) \]
In order to be able to perform anomaly detection, we must ensure that the sequence of incremental values \(\Delta_t\) can be approximated to a random walk. This can be verified by trying to fit an \(\text{AR}(1)\) process to the initial sequence of values \(y_t\) and checking whether the value \(|\alpha|\) is close to 1:
# Data lists for both delta-series
data_list_GDP <- list(
"y"=series_GDP,
"train_metrics_skip"=50,
"order"=1,
"N_train"=N_tot,
"N_test"=1 # No need to perform forecasting
)
data_list_CPI <- list(
"y"=series_CPI,
"train_metrics_skip"=50,
"order"=1,
"N_train"=N_tot,
"N_test"=1 # No need to perform forecasting
)
# Only need to inspect value of alpha for both series
to_save <- c(
"alpha"
)
# Fit AR(1) models
output_GDP <- fit_JAGS(
modelAR_string,
data_list_GDP,
to_save
)
output_CPI <- fit_JAGS(
modelAR_string,
data_list_CPI,
to_save
)
# Inspecting values of alpha
print(glue("Alpha value for GDP has posterior mean {output_GDP$mean$alpha} with standard deviation of {output_GDP$sd$alpha}"))
## Alpha value for GDP has posterior mean 0.859521726483674 with standard deviation of 0.030267447634924
print(glue("Alpha value for CPI has posterior mean {output_CPI$mean$alpha} with standard deviation of {output_CPI$sd$alpha}"))
## Alpha value for CPI has posterior mean 0.939175750757215 with standard deviation of 0.0200438333942664
Now that we know that the two series can be approximated reasonably well to random walks (\(\alpha_{GDP} = 0.86\) and \(\alpha_{CPI} = 0.938\)) we proceede with our anomaly detection.
Notice that the random variable \(\Delta_t\) has a fixed mean value \(\mu\), but a time dependant variance \(\sigma_t^2\), which we model as:
\[ \sigma_t^{-2} = \tau_t = \beta_1 + \beta_2 \delta_t \qquad \delta_t \sim Ber(p) \]
We can, therefore, define the following JAGS model-string which will fit the individual values of \(\delta_t\).
anomaly_detection_string <- "
model {
############################################################
######################## LIKELIHOOD ########################
############################################################
for(t in 1:N_tot) {
delta[t] ~ dnorm(mu, tau[t])
tau[t] <- beta[1] + beta[2] * d[t]
}
############################################################
########################## PRIOR ###########################
############################################################
# Beta coefficients, Bernoulli probability and mean
beta[1] ~ dexp(0.1)
beta[2] ~ dexp(0.1)
p ~ dbeta(1, 1)
mu ~ dnorm(0, 0.01)
# Delta value
for(t in 1:N_tot){
d[t] ~ dbern(p)
}
}"
# Obtaining series of differences
delta_GDP <- series_GDP[2:N_tot] - series_GDP[1:(N_tot-1)]
delta_CPI <- series_CPI[2:N_tot] - series_CPI[1:(N_tot-1)]
# Setting up parameters for JAGS
data_list_GDP <- list(
"delta"=delta_GDP,
"N_tot"=N_tot-1
)
data_list_CPI <- list(
"delta"=delta_CPI,
"N_tot"=N_tot-1
)
to_save <- c(
"d",
"mu",
"tau"
)
# Perform simulations
anomaly_output_GDP <- fit_JAGS(
anomaly_detection_string,
data_list_GDP,
to_save,
n_iter=10000
)
anomaly_output_CPI <- fit_JAGS(
anomaly_detection_string,
data_list_CPI,
to_save,
n_iter=10000
)
Once obtained both \(\{\delta_t\}_{t=1...N}\) sequences we can classify their values according to their expected values. Concretely, we classify a timestamp \(t\) such that \(\mathbb{E}[\delta_t] < 0.2\) as “high variance” (where strong oscillations may occur) and timestamps \(t\) with \(0.2 < \mathbb{E}[\delta_t] < 0.6\) as samples with “medium variance” (where we still might have an anomaly, but of weaker nature). All remaining samples can be considered stable and only subject to natural stochastic oscillations.
# Classifying anomalies for both series
anomalies_GDP <- classify_anomalies(anomaly_output_GDP, N_tot, title="GDP")
## Found 36 strong and 27 medium outliers out 305 total samples for GDP
anomalies_CPI <- classify_anomalies(anomaly_output_CPI, N_tot, title="CPI")
## Found 24 strong and 45 medium outliers out 305 total samples for CPI
# Displaying graphically the classification
delta_GDP <- anomaly_output_GDP$mean$d
delta_CPI <- anomaly_output_CPI$mean$d
display_classification(
x_axis,
delta_GDP,
anomalies_GDP$high_var_indices,
anomalies_GDP$med_var_indices,
anomalies_GDP$low_var_indices,
title="GDP"
)
display_classification(
x_axis,
delta_CPI,
anomalies_CPI$high_var_indices,
anomalies_CPI$med_var_indices,
anomalies_CPI$low_var_indices,
title="CPI"
)
Finally, we can plot both our time series displaying the points in time where high variances were observed: red vertical lines display timestamps where strong anomalies were detected, while orange dotted lines show medium anomalies.
# Zoom value: we only plot starting from this point on so as not to clog up the graph
start <- 150
visualise_anomalies(
x_axis,
series_GDP,
anomalies_GDP$high_var_indices,
anomalies_GDP$med_var_indices,
title="GDP",
start=start
)
visualise_anomalies(
x_axis,
series_CPI,
anomalies_CPI$high_var_indices,
anomalies_CPI$med_var_indices,
title="CPI",
start=start
)